﻿using System;
using System.Diagnostics;
using System.Text.RegularExpressions;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace Nature.Numerics
{
    [DebuggerDisplay("NxN[{Nva}][{Nva}]")]
    public class NxN : INxN
    {
        public class SimpleDecompositionFailedException : System.Exception { }

        #region Private Fields
        [DebuggerBrowsable(DebuggerBrowsableState.RootHidden)]
        private double[][] m_matrix;
        [DebuggerBrowsable(DebuggerBrowsableState.Never)]
        private int[] m_index;
        #endregion

        #region Constructors
        public NxN(int rank)
        {
            m_matrix = Extensions.AllocateMatrix<double>(rank, rank);
            m_index = new int[rank];
        }

        public NxN(double[][] matrix)
        {
            m_matrix = matrix;
            if (matrix == null)
            {
                throw new System.ArgumentNullException();
            }
            int size = matrix.Length;
            for (int i = 0; i < size; ++i)
            {
                if (matrix[i] == null)
                {
                    throw new System.ArgumentException();
                }
                if (matrix[i].Length != size)
                {
                    throw new System.ArgumentException();
                }
            }
        }

        public NxN(string text)
        {
            string pattern =
                @"(?xi) 
                    [+-]? \s*
                    (?:                    
                        (?: 
                            (?:\d+\.?\d*) | (?:\.\d+) 
                        )                            
                        (?:[edf][+-]?\d{1,3})?               
                    )                    
                ";
            double[] values = Regex.Matches(text, pattern)
                .Cast<Match>()
                .Select(m => double.Parse(m.Value))
                .ToArray();
            int rank = (int)Math.Sqrt(values.Length);
            if (rank * rank != values.Length)
            {
                throw new System.ArgumentException();
            }
            m_matrix = new double[rank][];
            for (int i = 0; i < rank; ++i)
            {
                m_matrix[i] = new double[rank];
                for (int j = 0; j < rank; ++j)
                {
                    m_matrix[i][j] = values[i * rank + j];
                }
            }
        }

        public static implicit operator NxN(string matrix)
        {
            return new NxN(matrix);
        }
        #endregion

        #region INxN Members

        public void Solve<T>(T rhs) where T : IVector
        {
            SimpleSolve(m_matrix, m_index, rhs);
        }

        #endregion

        #region IMatrix Members

        public double this[int i, int j]
        {
            get
            {
                return m_matrix[i][j];
            }
            set
            {
                m_matrix[i][j] = value;
            }
        }

        public void SwapRows(int i1, int i2)
        {
            throw new NotImplementedException();
        }

        public int Nva
        {
            get { return m_matrix.Length; }
        }

        public static implicit operator double[][](NxN matrix)
        {
            return (double[][])matrix.m_matrix.Clone();
        }
        #endregion

        public static void SimpleDecomposition(double[][] a, int[] index)
        {
            int rank = a.Length;
            for (int i = 0; i < rank; ++i) { index[i] = i; }
            for (int i = 0; i < rank - 1; ++i)
            {
                #region Pivoting
                {
                    int imax = i;
                    double scale = 0.0;
                    for (int ii = i; ii < rank; ++ii)
                    {
                        double[] rr = a[ii];
                        double current, maximum = 0.0;
                        for (int j = i; j < rank; ++j)
                        {
                            current = rr[j];
                            if (current < 0.0) { current = (-current); }
                            if (maximum < current) { maximum = current; }
                        }
                        if (maximum == 0.0)
                        {
                            throw new SimpleDecompositionFailedException();
                        }
                        current = rr[i] / maximum;
                        if (current < 0.0) { current = (-current); }
                        if (scale < current)
                        {
                            scale = current;
                            imax = ii;
                        }
                    }
                    if (imax != i)
                    {
                        index[i] = imax;
                        double[] r = a[i];
                        double[] rmax = a[imax];
                        for (int j = i; j < rank; ++j)
                        {
                            double tmp = r[j];
                            r[j] = rmax[j];
                            rmax[j] = tmp;
                        }
                    }
                }
                #endregion

                double[] mainRow = a[i];
                double mainElement = mainRow[i];
                for (int ii = i + 1; ii < rank; ++ii)
                {
                    double[] row = a[ii];
                    double factor = row[i] / mainElement;
                    row[i] = factor;
                    for (int j = i + 1; j < rank; ++j)
                    {
                        row[j] -= mainRow[j] * factor;
                    }
                    //if (row[ii] == 0.0)
                    //{
                    //    throw new SimpleDecompositionFailedException();
                    //}
                }
            }
        }

        public static void SimpleSolve<TVector>(double[][] a, int[] index, TVector rhs) where TVector : IVector
        {
            int rank = a.Length;
            for (int i = 0; i < rank - 1; ++i)
            {
                int imax = index[i];
                double tmp = rhs[i];
                rhs[i] = rhs[imax];
                rhs[imax] = tmp;
                for (int ii = i + 1; ii < rank; ++ii)
                {
                    double factor = a[ii][i];
                    rhs[ii] -= rhs[i] * factor;
                }
            }
            for (int i = rank - 1; i >= 0; --i)
            {
                double[] row = a[i];
                for (int j = i + 1; j < rank; j++)
                {
                    rhs[i] -= rhs[j] * row[j];
                }
                rhs[i] /= row[i];
            }
        }

        #region IDecomposible Members

        public void Decompose()
        {
            if (m_index == null)
            {
                m_index = new int[Nva];
            }
            SimpleDecomposition(m_matrix, m_index);
        }

        #endregion
    }

    //class Linearizible : NxN
}
